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ABSTRACT 

Linear stability calculations are presented for radially structured, vertically stratified, geometrically- 
thin disks with non-uniform entropy distribution in both directions. Polytropic equilibria are consid- 
ered but adiabatic perturbations assumed. The unperturbed disk has a localized radial density bump 
making it susceptible to the Rossby wave instability (RWI) . The linearized fluid equations are solved 
numerically as a partial differential equation eigenvalue problem. It is found that when the polytropic 
index is fixed and adiabatic index varied, non-uniform entropy has negligible effect on the RWI growth 
rate, but pressure and density perturbation magnitudes near a pressure enhancement increases away 
from the midplane. The associated meridional flow is also qualitatively changed from homentropic 
calculations. Meridional vortical motion is identified in the nonhomentropic linear solution, as well 
as in a nonlinear global hydrodynamic simulation of the RWI in an initially isothermal disk evolved 
adiabatically. Numerical results suggest buoyancy forces play an important role in the internal flow 
of Rossby vortices. 



1. INTRODUCTION 

Understanding the stability and evolution of radially 
structured disks is important for several astrophysical ap- 
plications. Protoplane tary disks are likely to have com - 
plex radial structure (|Terqueml 120081 : lArmitagel 120111 ) . 
Examples include the radial boundary between magneti- 
cally active and inactive regions of the disk ('dead zones', 
IGammid I1996D. and edges of gaps induced by a giant 



of 



planet (|Lin fc P apaloizoul ll986|) 

Local variations in the disk profile, which both 
the above examples involve, are vulner able to the so- 
called Rossby wave inst ability (RWI, lLovelace et al.l 
[1991 ILi et alJl2000L l200l . These studies assume a two- 
dimensional (2D) disk. The RWI is a linear shear insta- 
bility associated with an extremum in the potential vor- 
ticity profile of the disk, or a generalization thereof, and 
leads to local vortex formation in the nonlinear regime. 
This has been verified for b oth dead zone boundary 
and planetary gaps in 2D (e.g. Varmere fc Taggerll2006l; 
Lyra et all 120081 120091: ILi et all 120091: iLin fe Papaloizoul 
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These studies show that Rossby vortices are able to 
concentrate dus t particles, potentially assisting planetes- 
imal formation (Barge fe Sommerial ll995: Ina ba fc Bargd 
2006), which is of course crucial for planet formation. 
They can also interact strongly with planets, lead- 
ing to non-monotonic o rbital migration (lYu et alJl2010t 
ILin fc Papaloizoul l2"010). Although protoplanetary disks 
are thin, they are nevertheless three-dimensional (3D), 
so modeling these processes in 3D is necessary. 

Recently, the RWI has been demonstrated in 3D 
geometry in the c ontext of protoplanetary d isks 
(iMeheut et al.l l2010l l2012aHbllcl : (Umurhanl [20101: Err] 
I2012al|bl : iLvra fc Mac Low! 120121 ). These models have, 
however, employed a barotropic or nearly-barotropic 
equation of state. They can therefore be regarded 
as the thin-dis k version of the Papaloizou-Pringle in- 
stability (PPI, iPapaloizou fc Pringld H98l [19851 li~987l: 
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iGoldreich et~aT1 119861 : iNaravan et all 119871 ). originally 
discovered for 3D pressure-supported thick tori. It 
is clearly of interest to extend 3D RWI calculations 
to non-barotropic flow, which was one of the features 
that distinguished th e 2D RWI from the original PPI 
([Lovelace et al.lll999f ). 

Given that the RWI and PPI involve the same 
physics, that is, wave-coupling a cross co-rotation 
(|Goldreich et al.lll986HUmurhanll2010l ). it is worth point- 
ing out that the PPI has in fact been general ized to 
nonhomentropic tori. iFrank fc Robertson (|1988l) found 
that entropy gradients d id not significantly a ffect insta- 
bility growth rates, while iKoiima et all (|1989f ) concluded 
non-uniformity in entropy has similar effects as com- 
pressibility. They also found perturbations have weak 
vertical dependence, in agre ement with analytical argu- 
ments for homentropi c flow ([Papaloizou fc Pringfdll985l ; 
IGoldreich et~^[l986T) . 

In this work, we study what is essentially the nonho- 
mentropic PPI in rotationally-supported thin 3D disks, 
a geometry relevant to protoplanetary disks. This is 
equivalent to a n ext ension of the 2D RWI studies of 
lLovelace et all (|1999f ) to 3D, and we will adopt such 
nomenclature. 

We consider the problem in the linear regime. Al- 
though the role of the RWI in protoplanetary disks must 
be determined through nonlinear hydrodynamic simula- 
tions, linear calculations are nevertheless a useful way 
to study the instability at low computational cost. It 
is also important to have such calculations at hand for 
comparison with nonlinear simulations. 

Linear disturbances in 3D disks are governed b y com- 
plicated partial differential equations (lKata l2001l ). Even 
with a numerical approach, computing unstable modes is 
no simple task. One method is to evolve the linear equa- 
tions as an initial value problem ([Papaloizo u fc Pringld 
[19871: IFrank fc Robertson! [l988h and measure growth 
rates from data. For special disk equilibria, one can 
convert the p roblem to a set o f ordinary differen- 
tial equations ([Tanaka et al.l l2002t iZhang fc Lai 120061: 
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iMeheut et all l2012d; E3 |2012al hereafter ILH ). the 
derivation of which can be tedious. Thus, our study is 
also motivated by the desire to reduce this complexity 
when a numerical method is sought of. 

We pursue a numerical solution to the two-dimensional 
eigenval u e problem. This approach has been taken by 
IKoiimal (fl986l fl98l using finite-difference and finite- 
element methods. Inspired by previous studies, we 
employ finite differences in the radial direction and a 
seudo-spe ctral method to treat the vertical direction 
Linll2012q l. We formulate the linear problem with nu- 
merical implementation in mind, so that much of the 
algebra can be taken care of by the numerical scheme, 
should one choose to do so. 

This paper is organized as follows. In fj2]we list the gov- 
erning equations and describe the polytropic disk equi- 
libria under consideration. The linear problem is defined 
in $3] and the numerical method stated in $4l Linear 
simulations are presented in SJS] for disks with moderate 
values of the polytropic index. Disks with an isothermal 
background are considered in [J6] where a nonlinear hy- 
drodynamic simulation is also described. We summarize 
in ij7]with a discussion of important caveats and possible 
extensions to this study. 

2. DISK MODEL 

We consider a non-self-gravitating, inviscid fluid disk 
orbiting a central star of mass M* and adopt cylindrical 
co-ordinates (r, </>, z) centered on the star. The system is 
governed by the Euler equations: 



dp 

at 
dv 



v • (pv) = o, 
t . 



v ■ Vu = — Vp 

at p 



dt 



In s + v ■ V In s = 0, 



(1) 
(2) 

(3) 



where p is the mass density, v is the velocity field, p is the 
pressure and we refer to s = p/p 1 as the entropy, where 
the ratio of specific heats 7 is assumed constant. In the 
momentum equation, $* is the gravitational potential of 
the central star. Eq. [3] describes adiabatic evolution. 

A direct consequence of Eq. [I] and Eq. [2]is an equation 
for the vortensity (eVxo//), 



1 



— 7 = C • + — Vp x Vp, 



Dt 



(4) 



where D/Dt = dt + v ■ V is the Lagrangian derivative. 
The second term on the RHS is the baroclinic vorticity 
source. It is absent in barotropic flow for which p = 
p(p). In this work, we consider barotropic equilibria but 
generally non-barotropic disturbances, so the baroclinic 
term is effective in the perturbed state. 

2.1. Polytropic equilibrium 

The unperturbed disk is steady, axisymmetric and 
polytropic. That is 

p = Kp 1+1 ^, (5) 

where K is a constant and n is the polyt ropic index. 
We adopt the thin disk approximation (|L12D , so the den- 
sity field has the simple form p = p (r)(l — z 2 /H 2 ) n , 



where po(r) is the midplane density and H(r) is the disk 
thickness, po is specified indirectly by imposing a surface 
density profile £ cx r^ a B(r) where B(r) is a Gaussian 
bump at r = rp with amplitude A > 1 and width Ar 
(|Li et al.ll2000h . The aspect-ratio at r is parametrized 
as h = H(r )/r . 

The unperturbed velocity field is {v r ,v^,v z ) = 
(0, rO, 0) with O = f2(r) for barotropic equilibria and is 
given via centrifugal balance with gravity and pressure. 
Note that for a thin, non-self-gravitating disk the angu- 
lar velocity is nearly Keplerian, f2 ~ = J GM* / r 3 
where G is the gravitational constan t. 

The above setup is the same as in IL12I . and equations 
defining the equilibrium are listed therein. The limit n — > 
co corresponds to isothermal equilibria, and is treated as 
a special case in Sj6] 

Polytropic equilibria are ad opte d for simplicity and to 
allow direct comparison with IL12I . which considered ho- 
mentropic flow where L = 1 + 1/n = 7. Then Eq. [5]holds 
in the perturbed disk, replacing Eq. [3J Setting 7 / T 
gives a non homentropic dis k. 

Following Lov elace et al.l (|1999l) , it is convenient to de- 
fine the the following length-scales 



L, 



1 dlnp 

7 dr 

1 91ns 

7 dr 



H„ = 



1 dlnp 

7 dz 

1 d In s 

7 dz 



(6) 



(7) 



These are, respectively, the pressure and entropy length- 
scales in the radial and vertical directions, which depend 
on both r and z. Note that for polytropic equilibria, 
the entropy and pressure length-scales only differ by a 
constant multiplicative factor. 

2.2. Stability criteria 

We consider disk equilibria satisfying the Solberg- 
Hoiland criteria for stability against axisymmetric per- 
turbations: 



+ iV r 2 + N 2 Z > 0, k 2 N 2 > 0, 



(8) 



where k = r d(r ' fi )/dr is the square of the epicycle 
frequency and 



N 2 



LpL s 



N 



H P H S 



(9) 



are the radial and vertical buoyancy frequencies, respec- 
tively, and c s = ('jp/p) 1 / 2 is the adiabatic sound speed 
(|Tassou]|l200rl . We also define iV 2 = N 2 + N 2 . 

Our disk models satisfy the Rayleigh criterion k 2 > 0, 
which limits the surface density bump amplitudeQ- Then 
we require N 2 > 0, or stability against vertical con- 
vection, so the disk should be sub-adiabatically strati- 
fied with F < 7. For rotationally supported thin disks, 
\N 2 \ <C k 2 so the first Solberg-Hoiland condition is gener- 
ally s atisfied regardless of the equation of state (|Li et al.l 
2000). Note that N 2 increases with z, so we expect the 
disk to be more stable at larger heights. 

1 This also means that, by rescaling the density field, we can 
always make the Toomre stability parameter Qj- = c s k/ttGT, 3> 1, 
where c s is a typical sound-speed, to satisfy the assumption of a 
non-self-gravitating disk. 
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2.3. Instability criterion 

In t he original 2D RWI calculations, ILovelace et al.l 
(1999) found that when there is an extremum in the gen- 
eralized vortensity profile rj{r), where 



c 2 s 5p. The momentum equations give 



2ft£ 



n 



-2/ 72 



(10) 



the disk may be unstable to non-axisymmetric perturba- 
tions localized about the extremum. Here, IT = J_ pdz 
is the vertically integrated pressure and 72 is the adi- 
abatic index in the two-dimensional energy equation 
D{mi~^)/Dt = 0. 

To use E g. HUlin characterizing 3D disks, we use re- 
sults from IGoldreich et all (| 19861 ) to relate 72 and 7. 
iGoldreich et al.l studied linear disturbances in homen- 
tropic slender tori with a polytropic equation of state 
(Eq. [5]) . Assuming vertical hydrostatic equilibrium, they 
showed that the vertically integrated system has an ef- 
fective polytropic index of 712 = n + 1/2. If 72 = 1 + 1/712 
then 72 = (37 — l)/(7 + 1). This relation has bee n used 
by other authors (e.g . iLi et al.l 120001: |Klah3l200l. The 
assumptions made by IGoldreich et all do not strictly ap- 
ply to our case (nonhomcntropic equilibria and non-zero 
vertical motions) but their result will suffice for diagnos- 
tic purposes. 

The polytropic disk equilibria have II oc PqH and £ oc 
PoH, so the above definition gives 



77 oc 



K s (l-2r/ 72 )a-2(r-l)/ 7 2 

2fl 



For the adopted parameter values, a surface density 
bump corresponds to a local minimum in the general- 
ized vortensity, so that dn/dr ~ at r = tq. This is 
also close to a local min(«; 2 ). These minima act to 'trap' 
disturbances, leading to instability (|Li et al.H2000D . 



3. LINEAR PROBLEM 

We consider Eulerian perturbations to the above equi- 
librium in the form Re[<5p(r, z) expi(m</> + at)] and sim- 
ilarly for other fluid variables. Here, m is the az- 
imuthal wavenumber taken to be a positive integer and 
a = -co — \v is a complex frequency, where — oj is the real 
mode frequency and v is the growth rate. The co-rotation 
radius r c of a mode is such that mQ(r c ) — uj = 0, and the 
RWI is characterized by r c ~ r . For clarity, hereafter 
we omit writing out the time and azimuthal dependence 
explicitly. 

The goal is to obtain a partial differential equation 
(PDE) for the quantity W _ Sp/p. An explic it form of 
this equation is given by lKojima et al.l (|1989h . but our 
priority is the ease of solution implementation. By writ- 
ing individual equations in standard form — a sum of co- 
efficients multiplying differential operators — we can for- 
mulate the linear problem in convenient variables, then 
transform to the desired ones by redefining said coeffi- 
cients. These transformations can be done in the numer- 
ical code. 

We begin by writing down the linearized equations in 
terms of the intermediate variables W = pW and Q = 



pdv r ■ 
pSvj, 
p5v z ■■ 



i IdW 2mQ ~ x 
D \ dr r 



L P D 



Q, 



(ii) 



1 / k 2 dW ma. 



D \ 2VL dr 

i fdW Q 
a \ dz H v 



-W - 



2ttL p D 



Q, (12) 



(13) 



where a — a + mf2 is the shifted frequency and D = 
k 2 — a 2 . The linearized continuity equation is 

•- Q 1 3 , r \ 1171 r 9 , r \ r, AS 

io- -z H 7— (rpdvr) -I povf + — (pov z ) = 0, (14) 

c s r or r oz 

and the linearized energy equation is 

1 , . , 1 



ia(Q-W) =c 



(pSv r 



(pSv z 



(15) 



Inserting the momentum equations into the continuity 
and energy equations yield a pair of PDEs: 



W 




r dr \ L p D 



a dW 
L S D dr 



1 



L s LpD 



1 d 2 W 


2m 


d 


a dz 2 


T 


r 


dr 


1 d 


1 Q 


\ 




a dz 






dW 


2m 


n 


a 


dz 


rL s D 


V 2 

1 s 



r 2 D 



a H s Hp 



2mfl 
rLpD 



W 



Q = 0. 



= 0, 
(16) 

(17) 



Eq. 1161 — 1171 are the governing equations for linear distur- 
bances. 

Next, we transform to the co-ordinates (R, Z) — 
(r,z/H) so that the background disk structure is sep- 
arable. For example, the density field becomes p = 
po(R)g(Z). Then the unperturbed disk occupies a rect- 
angular domain since <?(±1) = 0. The governing equa- 
tions become 



d 2 W , d 2 W 

ai — g + 01- 



dR 2 



dZdR 



d 2 W , dW 
Cl -dZ^ + dl -dR 



dW 



dZ 



dW dW 



(18) 
(19) 



Explicit expressions for the coefficients are listed in Ap- 
pendix [A] We write the PDE pair for (W,Q) as 



^iW + yiQ = o, 

V 2 W + V 2 Q = 0. 



(20) 
(21) 



Since V2 is a multiplicative factor, we can eliminate Q 
between Eq. [201 — [2T1 to obtain an equation for W: 

(22) 



[Vi - Vi iy^ l v 2 )] w = vw = 0. 
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The operator V is obtained by updating the coefficients 
of V\ , so they have the same form. Finally, we substitute 
W = pW to obtain 



UW = 0. 



(23) 



Construction of V, and hence U , requires the evaluation 
of V\ (y^ 1 ]^) which involves radial and vertical deriva- 
tives of the coefficients in Eq. [19] In Appendix [B] we 
outline an alternative numerical approach which circum- 
vents the algebra. (This appendix also includes relevant 
formulae to redefine the PDE coefficients for the trans- 
formation V —> U.) 

The key dependent variable is W, but we also inter- 
pret results using Q = Q/p. We refer to W and Q as 
pressure and density perturbations, respectively. Then 
the entropy perturbation is naturally defined as 



S = W -Q. 



3.1. Boundary conditions 



(24) 



We consider disturbances radially confined about the 
density bump at r = ro, so the inn er and outer dis k 
boundaries play no significant role (|Umurhanl [2010). 
Hence, for simplicity we set OrW — at radial bound- 
aries. 

Pressure and density perturbations are assumed to be 
symmetric about the disk midplane. Henceforth we con- 
sider z > without loss of generality. The default upper 
disk boundary condition is vanishing Lagrangian pres- 
sure perturbation at Z = Z s , 



AP = SP + ^-\7P = 0, 



(25) 



where £ is the Lagrangian displacement (V refers to 
cylindrical co-ordinates). We call this the free boundary 
condition. The surface function Z s is assumed constant 
for simplicity. If Z s is the zero-pressure surface, then 
Eq. [25] can be satisfied automatically provided the per- 
turbations are regular there. In practice, though, we take 
Z s < 1 to avoid the disk surface (where entropy and its 
derivatives diverge, iZhuravlev fc Shakura 2007). Note 
that Eq. [25] together with Eq. OH imply TQ = -/W at 
the upper boundary. 
In some cases we adopt a solid upper boundary, 



5v z = Z s ^-5v r 
ar 



(26) 



meaning no flow perpendicular to the boundary (5v± = 
0). Upper disk boundary conditions are imposed explic- 
itly by replacing the governing equation with Eq. [25] or 
[2B|at Z = Z s . 

3.2. Baroclinity 

Before proceeding to solve the linear equations, it is 
useful to have a qualitative picture of the solution t o aid 
us in checking results. The main difference from IL12I 
is baroclinity. Here, we discuss expected effects of the 
baroclinic source term in Eq. 2] 

As we will often examine meridional flow, consider the 
azimuthal component of Eq. [4] which can source vor- 
tical motion in the (r, z) plane. When linearized, this 



baroclinic source term becomes 



(Vp x Vj>) 




where S = Q — jW/T and ' denotes differentiation with 
respect to the argument. We have utilized the barotropic 
background in obtaining Eq. 1271 In the discussion below, 
perturbations are regarded as real quantities. 

The RWI is characterized by non-axisymmetric pres- 
sure en hancements ra dially localized about the density 
bump (|Li et al.l 120011 ) . Assuming this is qualitatively 
unchanged in a nonhomentropic disk, let us denote the 
midplane co-ordinate of the center of one such enhance- 
ment as (ro,4>o). We will precisely define (j>o later. For 
now, consider the (r, z) plane at fixed <fi = </>o and about 
r = r . 

The distribution of S(R, Z) at the chosen azimuth can 
be anticipated as follows. First note that 



S 



w-s. 



For a pressure enhancement, (1 — 7 /T)W < 0. Recall 
s oc p r_7 so a density bump at ro corresponds to an en- 
tropy dip there. Then S(ro, 0) > 0, because the RWI has 
gathered material toward ro to increase the local pressure 
and density. The free boundary implies S(R,Z S ) = 0. 
Then it is reasonable to assume S(ro, Z) < 0. A similar 
argument can be made for the solid boundary. Finally, 
in moving radially away from ro, IS"! should decrease be- 
cause jthe RWI is radially localized. A simple picture is 
that S has a local minimum at (ro,0) and is negative 
or zero in this region. Then dzS > and 8rS > 
(d R S < 0) for R>r Q {R< r ). 

Where the disturbance is roughly two-dimensional 
(dz <C 1), which can be expected radially away from 
ro according to the above, the sign of Eq. [27] is dic- 
tated by that of (|). Even if dz = 0(1) in this re- 
gion, we expect Or ~ H _1 for a radially localized distur- 
bance. So the magnitude of (f) relative to (|) is of order 
\Hp' /Zp \(l — Z 2 ), which is small for the adopted disk 
models. Of course, this argument does not apply where 
dftS = 0, which occurs at Z = Z s and is expected close 
to r = ro. 

Thus, under the above assumptions we anticipate that 
away from the midplane but not close to the upper disk 
boundary, the baroclinic source term, Eq. 1271 is positive 
(negative) exterior (interior) to ro. Close to ro though, 
its sign the same as sgn(p ), which is typically negative, 
but not always, due to a density bump. 

4. NUMERICAL PROCEDURE 

The operator U can be written in the same form as V\ . 
A mat r ix repr esentation of such an operator is described 
in iLinl (|2012cD where details are given. We summarize 
here the main steps. 

The radial co-ordinate is discretized into Nr uniformly 
spaced grid points. Let Wi(Z) = W(Ri,Z) denote the 
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solution along the vertical line R — Ri. We set 

N z 

W i (Z) = ^w H il> h (Z/Z,) ) (28) 
fc=i 

where the basis functions ipk = ?2(fc-i) are even Cheby- 
shev polynomials of the first kind. Nz is the number 
of basis functions and the highest polynomial order is 

im« = 2(JV Z - 1). 

Radial derivatives in UW are replaced by by central 
finite differences, and we evaluate vertical derivatives 
exactly at the Nz non-negative Lobatto grid points of 
ri mn (Z/Z s ). This procedure performs the conversion 

UW = -> Uw = 0, (29) 

where U is a {NrNz) X (NrNz) block tridiagonal ma- 
trix and «; is a vector storing the NrNz pseudo-spectral 
coefficients Wki- 

The numerical problem is a set of linear homogeneous 
equations, U(a)w = 0. Non-trivial solutions exist if 
det (7 = 0. This is achieved by varying a using Newton- 
Raphson iteration. We only accept solutions where the 
reciprocal of the condition number of U is zero at ma- 
chine pre cision . The same method of solution was em- 
ployed in IL12I . 

4.1. Results visualization 

The pressure perturbation W is constructed from the 
pseudo-spectral coefficients Wki- We then calculate Q 
from Eq. [19] and velocity perturbations from Eq. [TT] 

ES 

We examine real perturbations about the vortex core 
(r,<f>) = (ro,</>o); where rn<fro = — axg[W(ro, 0)]. Setting 
4> — <f>o is equivalent to redefining a physical perturbation 

as 

X ->Re[X(r,z)W*(r ,0)], (30) 

where X represents W, Q, S or Sv, and * denotes com- 
plex conjugate. All perturbations are regarded as real 
hereafter. In practice (ro,<fio) is close to a local maxi- 
mum of pressure perturbation. The magnitude of X, as 
redefined above, is arbitrary but its sign is not. 

As an empirical measure of flow three-dimensionality, 
we compare vertical and horizontal motions near the 
bump radius using (0 m ), where 

° 2m = 6v?+svr (31) 

and (•) denotes averaging over R € [0.8, 1.2] ro and Z G 
[0, Z s ] at 4> = 4> . 

5. LINEAR SIMULATIONS 

We adopt units such that G = M» = 1. Our main cal- 
culations are summarized in Table [T] For these runs the 
computational domain is R £ [0.4, 1.6]r , Z e [0,Z S ] = 
[0, 0.9], and a = 0.5 for the power-law part of the surface 
density profile. The bump radius, amplitude and width 
are set to ro = 1, A = 1.4 and Ar = 0.05ro, respectively. 
We consider modes with m , = 3 unless otherwise stated. 
Slightly different setups are employed in fjS] to explore 
the isothermal limit. 

The new parameter for nonhomentropic disks, com- 
pared to homentropic flow in IL12I . is the adiabatic index 




o.o L , , i , , , i , , , i , , , i , , , i , , J 

0.4 0.6 0.8 1.0 1.2 1.4 1.6 

r/r a 

Fig. 1. — Equilibrium profile for a nonhomentropic disk with 
n = 1.5 and 7 = 2.5 (case 3). The generalized vortensity (top) and 
k 2 + N 2 at three heights (bottom) are shown. 



TABLE 1 

Summary of main linear simulations 



Case 


r 


7 


BC+ 




u/mQo 


io 2 is/n 


<M 








h -- 


= 


0.14 









1.67 


1.67 


AP = 





0.9941 


10.74 


0.33 


1 


1.67 


1.8 


AP = 





0.9937 


10.80 


0.36 


2 


1.67 


2.0 


AP = 





0.9931 


10.86 


0.39 


3a 


1.67 


2.5 


AP = 





0.9919 


10.99 


0.44 


3b 


1.67 


2.5 


8vj_ = 





0.9911 


11.34 


0.41 


4 


1.67 


3.0 


AP = 





0.9910 


11.07 


0.47 








h 




0.2 






5 


1.4 


1.4 


AP = 





0.9923 


16.66 


0.24 


6 


1.33 


1.4 


AP = 





0.9917 


13.81 


0.31 


7 


1.29 


1.4 


AP = 





0.9912 


11.38 


0.34 


8 


1.25 


1.4 


AP = 





0.9909 


9.246 


0.36 



t Boundary condition at Z = Z s . 



7. We therefore focus on examining the effect of entropy 
gradients due to 7 ^ T. Cases — 4 have fixed polytropic 
index n = 1.5, and therefore identical background den- 
sity and velocity profiles, but variable adiabatic index 
7 > 5/3. Cases 5 — 8 have fixed adiabatic index 7 = 1.4, 
but variable polytropic index n > 2.5. 

An example of nonhomentropic equilibrium, with n = 
1.5 (r = 5/3) and 7 = 2.5, is shown in Fig. Q] (case 
3). The generalized vortensity and k 2 + N 2 are plotted. 
As expected for a density bump, the generalized vorten- 
sity has a local minimum at r = ro- It corresponds to 
min(K 2 /f^) = 0.43. The increase in k 2 + N 2 with re- 
spect to height is due to N 2 (since N r ~ hN z near the 
upper boundary). Note that N 2 > fl 2 for \z\ > 0.7 H in 
this case. 

The discretized problem is solved with standard ma- 
trix routines provided in the LAPACK package. The de- 
fault resolution is (Nr,Nz) — (512,12), corresponding 

tO ^max — 22. 
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Fig. 2. — Real perturbations for the homentropic case (T = 
7 = 1.67). The pressure perturbation W at three heights (top) and 
meridional velocity perturbation (bottom) near the vortex core are 
shown. 

5.1. Homentropic reference case 

For comparison purpos es, we reproduce the fiducial ho- 
mentropic calculation in IL12I by setting 7 = 5/3 (case 
0). Then W = Q since L~ l = H' 1 = (Eq. [15}. 
This also serves as a test for our numerical method. 
The eigenfrequency and pert urba tions shown in Table 
Q] and Fig. [2] agrees well with IL12I . In co-rotation region 
R G [0.8, 1.2] ro, W is nearly independent of height and 
the vortex core has upwards motion. 

5.2. Nonhomentropic example 

We now examine case 3a with 7 = 2.5, T — 1.67. The 
cigenfrequecy a is close to case 0, but the growth rate is 
slightly larger in the nonhomentropic disk. 

Fig. [3] shows the the pressure, density and entropy 
perturbations at several heights. Near ro, pressure and 
density perturbations increase with height, unlike the ho- 
mentropic case where W has weak z-dependence. Non- 
homentropic disks generally have W ^ Q, as shown in 
Fig. [3] The difference between W and Q at the midplane 
is due to background radial entropy gradients Lj 1 since 
H- 1 (r,0) = 0. 

At co-rotation, Q increases with height faster than the 
pressure perturbation, which results in a negative en- 
tropy perturbation. This is consistent with the require- 
ment S = (1 — 7/r)T / F at the upper disk boundary. It 
is clear that S has a stronger vertical dependence than 
either W or S. 

We might have expected the above result on physi- 
cal grounds. The homentropic case indicate upward mo- 
tion at the vortex core. If a positive (stable) vertical 
entropy gradient is introduced, then a fluid element dis- 
placed upwards should increase its density compared to 
the surrounding background, i.e. Q > 0, and this should 
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Fig. 3. — Pressure (top, W), density (middle, Q) and entropy 
(bottom, S) perturbations for the nonhomentropic case 3a (r = 
1.67, 7 = 2.5). 

become more positive with height because vertical veloci- 
ties increase in magnitude with height. The pressure per- 
turbation is not expected to change as rapidly, because 
the fluid element can establish pressure equilibrium with 
its surroundings. 

5.2.1. Entropy perturbation 

We plot the entropy perturbation S at z — and z — 
0.8H in Fig. @] The figures are overlaid by the perturbed 
horizontal flow, which are similar at both heights. The 
anti-cyclonic flow pattern i s comm only found in previous 
studies (e.g. ILi et all 120001 [2001 . Entropy gradients of 
this magnitude do not affect this characteristic feature of 
the RWI. 

Therefore, we could have inferred some of the features 
in Fig. 2] without solving the fluid equations. Consider 
the linearized energy equation near co-rotation where 
fx — — iz/ (which is close to ro), 

Ss v^Sv-Vs, (32) 

and v > for a growing mode. Eq. [32] is only valid 
within a small distance e <C ^/Im-figl from ro. In this 
example, J//|mf2 | ~ 0.02ro. 

The midplane entropy increases globally in the radial 
direction with a dip at ro, so d r s > there. This natu- 
rally implies that inward (outward) radial flow for <j) < 4>q 
((f> > <po), i.e. anti-cyclonic motion, brings about a local 
entropy increase (decrease) near r at z = 0. We recog- 
nize the qualitative similarity between the midplane flow 
pattern in Fig. [4] and horseshoe turns induced by an em- 
bedded planet. Entrop y advection then leads to l arge ra- 
dial entropy gradients (Paardckooper et al. 2010), which 
can be seen in Fig. on cither side of r at the mid- 
plane. This gradient is, of course, growing exponentially 
in time, so it may be important even within the linear 
regime. 

In the vertical dimension, if we assume the flow at 
(ro, 4>a) is unchanged from the homentropic case (i.e. up- 
ward) , then since the background entropy increases with 
height, the local Eulerian entropy perturbation at the 
vortex core must become negative away from the mid- 
plane, as observed. 
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Fig. 4. — Entropy perturbation in midplane (left) and near the 
upper disk boundary (right) for the nonhomentropic case 3a. Ar- 
rows show the perturbed velocity field projected onto this plane. 

_ Related to the entropy perturbation is the quantity 
S = Q — "fW/T. Its distribution shown in Fig. [5] agrees 
with expectations made in §3.2| namely it is mostly neg- 
ative, with a local minimum at the vortex core. 




Fig. 5. — Map of the quantity S = Q — fW/F for case 3a (the 
real perturb ation at <p = 4>0 is shown). S appears in the baroclinic 
ter m in Eq. 1271 as well as the expression for vertical velocity in 
Eq. 1331 The local minimum ne ar (r p , 0) can be expected without 
solving the linear problem (see i|3.2ll . 



5.2.2. Meridional vortical flow 

Fig. [5] shows the perturbed velocity field in the (r, z) 
plane, with a map of the baroclinic source term defined 
in i j3.2l The flow pattern is similar to the homentropic 
case in that it is still converging toward r*o, and verti- 
cal motion is predominantly upwards there. However, 
there is a notable difference from the homentropic case 
— vortical motion (of positive azimuthal vorticity) cen- 
tered about (r, z) = (1.02rn, 0.5H). It coincides with 
a region where the azimuthal baroclinic source term is 
positive. Note that the sign of the baroclinic source is 
roughly consistent with expectations in 




Fig. 6. — Meridional flow in the nonhomentropic case 3a. The 
azimuth taken for this slice is <j> = (f>o- The cont ours show the 
baroclinic source term for azimuthal vortensity (Eq. 1271 multiplied 
by p). The arrows show the perturbed velocity field projected onto 
this plane. 




Fig. 7. — Pressure (W , top) and density (Q, bottom) pertur- 
bation for the m = 5 mode in the nonhomentropic case 3a. The 
meridional flow is also shown. 

5.2.3. m = 5 

The meridional flow varies with m. Fig. [7] shows the 
7?i = 5 solution for the setup of case 3a. We focus on 
the region R £ [0.9, l.l]r bec ause higher-m modes ar e 
not as well- localized as low-m (jLin fc PaDalo izoull2011| ). 
It displays stronger vortical motion than the fiducial run 
with m = 3, even though the growth rates are similar 
(v/mVtQ = 0.1051 for m = 5). The pressure and density 
perturbations have noticeable vertical structure, with W 
typically increasing away from the midplane. This qual- 
itatively differs from homentropic cases. 

5.3. Solid upper boundary 

In the above example, it is perhaps not surprising that 
entropy perturbations became more negative away from 
the midplane, because the free boundary condition de- 
mands \Q\ > \W\ at Z = Z s . 

We have re-calculated this mode with a solid upper 
disk boundary (case 3b). Numerically, this condition 
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forces W ~ Q at Z = Z s . Fig. [5] shows the ratio of 
pressure to density perturbation. The entropy pertur- 
bation at intermediate heights is still typically negative, 
suggesting this to be an intrinsic feature of the instability 
in these disk models. The flow pattern is very similar to 
case 3a. 
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Fig. 8. — Ratio of pressure to density perturbations at </>q, for the 
nonhomentropic disk with solid upper disk boundaries (case 3b). 
A ratio above unity implies positive entropy perturbation. 



5.4. Effect of 7 on vertical flow 

When 7 ^ T, the presence of buoyancy forces is ex- 
pected to modify the vertical flow associated with the 
RWI. In Fig. |9] we compare the vertical velocity at the 
vortex core for a range of 7. 




m 0.0 



Fig. 9. — Normalized vertical velocities at the vortex core (ro, (f>o) 
as function of 2, for several values of 7 with fixed V = 1.67. The 
dash-dot line employed a solid upper disk boundary, other cases 
use the free boundary condition. 

As we increase 7, the magnitude of vertical flow in- 
creases, with an increasingly complicated z-dependence. 
In the homentropic case (7 = 1.67), Sv z is essentially 
linear in z, con sistent with the analytical models of 
lUmurhanl (|2010f ). For 7 > T, near z — it is still linear 
in z, but away from the midplane the increase in 8v z lev- 
els off. Sv z increases again at large z, but this is likely 
a boundary effect as indicated by a comparison between 
the dashed and dash-dot lines which employ different up- 
per disk boundary conditions. 

At (ro, <po), the vertical velocity is roughly 

1 \dW 2nZT 
~Vtt [W + 7(1 - Z 2 ) 
t ^ 



Sv z 



(q 



W 



(33) 



The first term (f) represent pressure forces and is present 
for all values of 7. For homentropic flow, (f ) is the only 
source of vertical motion, and in this case W decreases 
with height. The second term (J) is only present if 7 7^ T. 
Recall the quantity Q — jW/T = S defined in H3.21 where 
it appeared as a baroclinic source term and we argued 
S < at the vortex core (see also Fig. [5]). Then at the 
vortex core, (|) contributes positively to 5v z along the 
vertical direction, but vanishes at endpoints. 

In the nonhomentropic example (case 3a, 7 = 2.5) the 
function W increases with height at the vortex core (Fig. 
[3]), implying (f) contributes negatively to Sv z . The fact 
that we observe positive vertical velocity shows that (J) 
is typically larger in magnitude than (f). 

5.4.1. The role of 

Notice even when 7 is only slightly larger than F, the 
vortex core vertical velocity is quite different from the ho- 
mentropic case (i.e. case 1 with 7/r = 1.08 in Fig. [5]). To 
see the role of entropy gradients, or equivale ntly th e effect 
of non-zero buoyancy frequency, we follow iKatol (pOOlh 
and make the following approximations. For generality, 
below we shall not specialize to a polytropic background. 

Consider a height at which H^ 1 ^> L^ 1 , which is gen- 
erally true away from the midplane of a thin disk. Fur- 
thermore, suppose radial velocities are not much larger 
than vertical velocities in the region of interest (co- 
rotation). Then we can neglect the Sv r term in the lin- 
earized energy equation, and eliminate Q between Eq. 
[131 and Eq. [IS] to obtain, 



6v~ 



dW 

dz 



dlnp 

dz 



W 



5v z , (34) 



which is lKatof s Eq. 21 evaluated at co-rotation. Because 
v <C r^o for the modes considered and N z ~ Q away from 
the midplane, for nonhomentropic flow we have N 2 jv 2 ^> 
1 and should expect the balance 



Sv z 



dW 



N 2 dz 



v 

'N 2 



<91np 

dz 



H, 



W, N^O 



w 



(35) 

near co-rotation radius. Comparing the coefficient of 
d z W in this expression to that in the strictly homen- 
tropic case, where 

1 dW 

v dz 



Sv z 



N? = 0, 



we see d z W is less significant when N 2 ^ 0. Away from 
the midplane of a nonhomentropic disk, d z W becomes 
unimportant compared to the second term on the RHS 
of Eq. [35j which is just buoyancy (and does not explic- 
itly depend on 7). So the origin of vertical motion is 
qualitatively different between homentropic and nonho- 
mentropic flow away from the midplane, as suggested by 
numerical results in the previous section. 



5.5. 



ible r 



i 



Fixed 7, ■ 

We now fix the adiabatic index to 7 = 1.4, as is typ- 
ical for accretion disk models. Then we require n > 2.5 
for axisymmetric stability. With other parameters fixed, 
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increasing n would decrea se th e bump in disk thickness 
and reduce growth rates ()L12| ). To avoid potential nu- 
merical issues associated with small \a\ at co- rotation, 
we adopt h — 0.2 for cases 5 — 8, so that growth rates 
remain O(0.1il ). 

Table Q] shows that by setting 7 ^ T = 1.33 (case 
6), thereby introducing entropy gradients, (0 m ) has in- 
creased from the homentropic case 5. This is consistent 
w ith t he trend in cases — 4. 

IL12I found that when n is increased but other parame- 
ters fixed, the flow at the vortex core became less three- 
dimensional. For cases 5 — 8, we find the average value of 
9 m , when taken over R £ [0.98, 1.02]r , is 0.46, 0.63, 0.61 
and 0.56 for n — 2.5, 3.0, 3.5 and 4.0, respectively. The 
flow at (ro,</>o) m fact becomes more three-dimensional 
when it is nonhomentropic although n has increased (case 
5 — > case 6). 

The small decrease in the above values of three- 
dimensionality at (ro,</>o) from n = 3.0 to n = 4.0 is 
likely related to increased radial flow across r associ- 
ated with vortical motion in the (r, z) plane. Cases 6 — 8 
display similar dependence of 5v z on z as the nonhome- 
ntropic example (case 3a, see Fig. [9]). 

6. ISOTHERMAL LIMIT 

We examine the limit where the unperturbed disk is 
isothermal, but with adiabatic index 7 = 1.4. These cal- 
culations provide another check on our numerical results. 

6.1. Large polytropic index 

We first consider setting n = 10 to produce an almost 
radially isothermal equilibrium with p oc p . This al- 
lows us to use the numerical code as set up for polytropic 
equilibria without modification. We also adopt A = 2.5 
and h = 0.25 for reasons given in £15.51 The relatively 
large aspect-ratio does not violate the thin-disk approxi- 
mation as large n implies the density decays rapidly away 
from the midplane. Also because of this, we set the upper 
disk boundary at Z s = 0.6 to avoid very low densities. 

For this setup we obtained uj/mVlQ = 0.9883, v/VLq = 
0.1375 and (6 m ) = 0.35. The top panel of Fig. EUshows 
the meridional flow at the vortex core. The vortical mo- 
tion is distinct and more apparent than case 3a, despite 
the smaller value of 7/r in the present case. However, 
apart from this difference, the solution is qualitatively 
similar to case 3a. 

6.2. Strictly isothermal equilibrium 
Modifications to our standard setup are required to 



treat disk equilibria with p 
the constant sound speed c ls 



c loP ( r = !): where 

ft-iso^o( r Ao) 3 ^ 2 is the isothermal scale-height, and hi BO is 
the characteristic aspect-ratio at r$. The dimensionless 
vertical co-ordinate is now Z — z/H- lso . The isothermal 
atmosphere is exponential, g{Z) = cxp (-Z 2 /2), so there 
is no surface, but in practice we choose a finite vertical 
domain. 

In the linear code we simply replace expressions for the 
entropy and pressure length-scales by those correspond- 
ing to the isothermal disk, the function H by iJj so and 
g(Z) becomes the Gaussian above. We set parameter 
values to mimic the large-n polytrope in the previous 
section: Z s = 3 and hi BO = 0.05, so that for both setups 




Fig. 10. — Perturbed meridional flow at (j> = (fig for a n = 10 poly- 
tropic disk equilibrium (top) and a strictly isothermal equilibrium 
(bottom). 




Fig. 11. — Pressure (top, W), density (middle, Q) and entropy 
(bottom, S) for a globally isothermal background. 

the pressure is reduced by roughly the same factor in go- 
ing from the midplane to the upper boundary, and the 
midplane temperature is roughly equal at r$. 

We obtain w/mO = 0.9860, v/Q, = 0.1008 and 
(6m) — 0.39. The perturbations plotted in Fig. [TT] are 
similar to case 3a, so we expect these are features of 
the RWI in nonhomentropic flow, rather than associated 
with particular parameter values. The perturbed merid- 
ional flow shown in Fig. [10] (bottom panel) is in quali- 
tative agreement with the large-n polytrope. The result 
is, however, qui te different to i s otherm al linear perturba- 
tions, for which iMeheut et al.l (|2012d ) found the vertical 
velocity appears to have a node at ro (i.e. the fluid col- 
umn is hydrostatic there). Note that both 7/r and the 
growth rate are slightly smaller than the nonhomentropic 
case 3a, but here the vortical motion is more prominent. 
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6.3. A nonlinear simulation 

We have also performed global 3D hydrodynamic 
simulations using the ZEUS-MP finite-difference code 
(|Haves et a l. 2006). As the focus of this work is the linear 
problem, though, we defer a full discussion of these non- 
linear simulations to a follow-up paper. Our priority here 
is to verify the vortical motion in the meridional plane, 
which appears characteristic in the linear RWI solution 
for nonhomentropic flow. 

6.3.1. Setup 

We use spherical polar co-ordinates (r sp h,0, <f>) to de- 
scribe the disk, taken to be initially isothermal as de- 
scribed above. The computational domain is r sp h £ 
[0.2,2.0]r , 9 e [0 min , tt/2], <p £ [0, 2vr] and is divided 
into (512,48,512) zones, with tan (n/2 — 8 min ) = 3h iso 
and ro = 10. The grid is logarithmically spaced in ra- 
dius and uniformly spaced in the angular co-ordinates. 
Boundary conditions are outflow in r sp h, reflection in 9 
and periodic in <j). Additional damping to meridional ve- 
locities nea r radial boundaries are employed to reduce 
reflections (|de Val-Borro et al.ll2007l) . 

After some experimentation, we found it was most 
convenient to start with a smooth disk. In this case, 
a surface density £ oc r~ 3 / 2 , and tape red toward the 
inner boundary (as used in iLinl [2012bO . This reduces 
numerical transients associated with initialization with 
large radial gradients. We introduce the density bump at 
r = ro via source terms in the mass, momentum and ther- 
mal energy equations, over a time-scale of 10Po, where 
Pq = 2n/ilk(ro)- At t = 10Po we add small-amplitude 
random radial velocity perturbations. 

We choose the bump amplitude A = 1.2 5 and isother- 
mal asp ect-ratio hi SO — 0.1, as employed bv lMeheut et al.l 
(2012cJ) so that we can check our results against theirs. 
We measure perturbations with respect to azimuthally 
averaged hydrodynamic quantities at t = IOPq. 

6.3.2. Results and comparison to linear flow 

We focus on the earliest stage of the instability, when 
perturbation amplitudes are small so comparison with 
linear calculations can be made. Fig. [T2] shows the 
snapshot to be examined. A m = 4 mode has devel- 
oped. Notice the double-peak in density perturbation, 
which is also present in Fig. [TT] Using the method 
described in Appendix [Cj we estimated the m = 4 
mode growth rate and frequency to be z //flp ~ 0.189 
and u/ mQg ~ 0.990, in agreement with iMeheut et al.l 
(|2012cf ). Although they assumed barotropic perturba- 
tions, whereas we simulate adiabatic evolution, our linear 
calculations indicate growth rates are largely unaffected 
by entropy gradients (Table [TJ . 

We have also computed this mode using the linear code 
(as modified for strictly isothermal equilibria, with solid 
upper boundary). We obtain growth rate and mode 
frequency z^/^o = 0.1937 and w/mf2 = 0.9896, re- 
spectively. This is close to the nonlinear simulation. 
IMeheut et al.l (|2012cD suggested smaller growth rate in 
the latter may be due to numerical viscosity. Fig. 1131 
compares the density perturbation Q computed from 
the hydrodynamic simulation and linear code. They are 
broadly consistent. The linear code also produces a bias 
toward the over-density ahead of the vortex core at the 
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Fig. 12. — Nonlinear hydrodynamic simulation of the RWI in a 
nonhomentropic 3D disk, initially isothermal but evolved adiabat- 
ically. The axes are in units of ro. The relative density pertur- 
bation near the midplane, scaled by 100, is shown. This quantity 
is proportional to the Q used in linear calculations. The snapshot 
corresponds to the onset o f in s tab ility. The line drawn defines the 
vortex azimuth <j>o in Fig. 1131 — II II 

midplane. Away from the midplane, the center of the 
anti-cyclonic motion has shifted downstream. This shows 
that, even within the linear regime, the vortex is not 
columnar and has non-negligible vertical structure in the 
density perturbation (by comparing the two heights in 
Fig. [TJ. 

We compare meridional flows in Fig. 1141 The nonlin- 
ear simulation also produce vortical motion in the same 
sense as the linear calculation. The asymmetry of the 
pressure perturbation about ro is captured by the linear 
code as well. Disagreement toward the upper boundary is 
not unexpected, since the linear code assumes the upper 
boundary is at a constant number of scale-heights above 
the midplane, whereas the spherical grid imposes con- 
stant opening angle, and the RWI is a global instability 
in the vertical direction. However, both plots indicate W 
increases away from the midplane in the region exterior 
to r . 

7. SUMMARY AND DISCUSSION 

In this paper, we have examined the linear stability 
of radially structured three-dimensional disks with non- 
uniform entropy distribution. These calculations may be 
consid ered as an exte nsion to the 2D Rossby wave insta- 
bility (|Li et al.l I2000T ) by adding the vert ical dimension, 
or to the barotropic RWI calculations of IL12I by adding 
an energy equation with a simpler numerical method. 

We adopted polytropic disk equilibria so that the 
magnitude of entropy gradients can be conveniently 
parametrized by A7 = 7/T, and we focused on the effect 
of A7 > 1. When the background density and veloc- 
ity field is fixed through T, we found increasing A7 has 
negligible effect on the instability growth rate. However, 
pressure and density perturbations increase with height, 
and the meridional flow associated with the vortex core 
is qualitatively changed, with the introduction of vortical 
motion. We also found that the vertical velocity at the 
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Fig. 13. — Normalized density perturbation, Q, associated with 
the RWI computed from a nonlinear hydrodynamic simulation 
(left) and the linear code (right), at the midplane (top) and at 
2 scale- heights away from the midplane (bottom). The perturbed 
velocity field is also shown. The azimuthal wavcnumber is m = 4. 

vortex core is no longer linear in z, as for homentropic 
flow. 

In our second set of experiments, we fixed 7 and de- 
creased r. We found that by making the flow nonho- 
mentropic, the co-rotation region became more three- 
dimensional, despit e the decrease in growth rate. This 
result is opposite to lL12l where lowering T made the flow 
less three-dimensional. This implies that entropy gradi- 
ents play an important role in the vertical structure of 
the perturbations. 

We also considered isothermal equilibria. A linear cal- 
culation with r = 1.1 and one with a strictly isothermal 
setup (r = 1) were consistent. Both produced promi- 
nent meridional vortical motion. In order to verify this 
feature, we ran a nonlinear simulation of the RWI in 
an initially isothermal disk, but evolved adiabatically. 
We indeed identified said vortical motion. Keeping in 
mind that the setup for linear and nonlinear simula- 
tions were not identical (e.g. numerical grid, boundary 




Fig. 14. — The perturbed velocity field projected onto the merid- 
ional plane at the vortex azimuth <f>o , associated with the RWI cal- 
culated from a nonlinear hydrodynamic simulation (top) and the 
linear code (bottom) . A map of the normalized pressure perturba- 
tion is also shown. 

treatment), similarities between them, such as mode fre- 
quency, growth rate and horizontal flow, are satisfactory. 

Vortical motion in the meridional plane thus appears 
characteristic of the linear RWI in nonhomcntropic disks. 
Whether or not this is significant for the vortex evolu- 
tion can only be answered by detailed long term nonlin- 
ear simulations. If this vortical motion is present in the 
nonlinear regime then it may prevent dust particles from 
reach ing the disk surface, which occurs for homentropic 
flow (|Meheut et alJl2012bf ). 

However, given this meridional vortical motion is ab- 
sent in the homentropic linear solution, it may eventu- 
ally vanish because of entropy mixing, if no mechanism 
is present to maintain entropy gradients. For example, 
the background entropy increases with height but the 
linear entropy perturbation, which is growing exponen- 
tially, becomes more negativ e with height. On the other 
hand, iMeheut et al.l (|2012d ) observed strong meridional 
vortical motion in their homentropic hydrodynamic sim- 
ulations. We conclude they are of nonlinear origin. 

In the linear solutions, we often observe perturbation 
magnitudes increase away from the midplane in nonho- 
mentropic diskf0 (e.g. Fig. Q3]). Then the RWI may 
not be as robust against vertical boundary conditions 
as it is to radial boundary conditions. This could pose 
difficulty for the RWI to develop in dead zones of real 
protoplanetary disks, which are expected to be confined 
from above and below b y magnetically turbulent layers 
(|Oishi fc Mac Lowll2009D . The vertical boundary condi- 
tion set by these layers may or may not be compatible 

2 This reminds us of t he off-midplane vortices discovered by 
Barranco & Marcus (2005) in nonlinear local simulations, but the 
setup considered in that study is very different from the present 
work. Nevertheless, in both cases the vertical entropy gradient is 
stabilizing away from the midplane. 
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with the linear RWI solution. 

7.1. Caveats and outlooks 

One trade-off for the simplicity of our numerical 
method for linear simulations is that a trial eigenfre- 
quency must be guessed. This is not a significant ob- 
stacle for the problem at hand, b ecause previou s RWI 
studies provide an important guide t|Li et al. 2000). Oth- 
erwise, zeros of the complex function V(a) = dct U need 
to be located with more rigoro us methods (e.g. iKoiimal 
119861: Ide Val-Borro et all [20071) . We have also exploited 
previous findings that the RWI is predominantly two- 
dimensional, which enabled the use of a small number 
of basis functions. However, there could exist parameter 
regimes where the RWI has significant vertical structure, 
rendering our solution method inefficient. 

Our conclusions are limited to polytropic backgrounds. 
While this was convenient for numerical experiments, it 
is an over-simplification of protoplanetary disks, which 
are expected t o have complicated vertical structure 
(|Terqueml I2008D . In particular, we found that entropy 
gradients plays a role in the vertical structure of the RWI, 
and even a small entropy gradient can noticeably mod- 
ify the vertical flow ( §5.4. 1|) . Thus, a realistic model for 
entropy evolution is needed. 

It would also be of interest to generalize the calcula- 
tions to baroclinic equilibria, for which d z Q ^ 0. This 
may well be the case when the equilibrium pressure de- 
pends on both the density and temperature. Compli- 
cations from baroclinic instabilities may arise, however 
(jKnobloch fc Spruit! [1981 lUmurhad 12012 ). 

We have neglected gas self-gravity in this study. Our 
models therefore assume that the Toomre parameter is 
much larger than unity in both the unperturbed and per- 
turbed states. Previous studies have found higher m 
RWI modes are favored when disk self-grav ity is included 
(|Lvra et al.1 [20(51 iLin fc Papalorzoul l20lTh . Recent 3D 
simulations of the RWI in a locally isothermal disk show 
that vertical self-gravity can noticeably enhance the den- 
sity perturbation near the midplan e, even tho ugh the 
initial disk was considered low mass (|Linll2012bD . 

In principle, one can express the Poisson integral as a 
matrix operator and incorporate it into our formalism. 
This problem is further complicated by the need of a nu- 
merical solution to the equilibrium equation s describing a 
radially structured, self- gravitating 3D disk (|Mutoll20Tll ). 
Such a linear calculation is beyond the scope of this pa- 
per, but will be inevitable for understanding the RWI 
in 3D self-gravitating disks. Perhaps a simpler start- 
ing point, to gain first insight, is direct hydrodynamic 
simulations. This is indeed the approach taken in our 
follow-up paper. 



3 In fact, ba r oclini c tori were briefly considered by 
IFrank fc Robertson! (fT988T >. 
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APPENDIX 
PDE COEFFICIENTS 

In (R, Z) co-ordinates, the coefficients for the PDE pair (Eq. IT81 — FT9]) with dependent variables (W, Q) are 



ai = 1, bx = -2Z- 



'H' 



.w_ 



H 



di 
fi 



LpRo 



ei = 
D 

' ^2 



H 

Z H' 
T p ~H 



c\ = Z 



ml* 

D 
D 



D 



di = 



In 




a 2 HH p ' 



D dH- 
d 2 !! dZ 



dL- 1 
dR 



2m,n 
Rd 



In — 



m 



In — 



and 



h 



L S D' 
1 



e2 = 
1 



DL s L p 



Za H' 
~Lj5~H 
1 



1 



h = 



H' dL~ l 
Z p_ 

H dZ 
2mfl c 



(Al) 



L S RD 



aH p H s 



(A2) 



Note that these coefficients are expressed in terms of pressure, entropy length-scales and the adiabatic sound speed. 
Although H has the physical meaning of the polytropic disk thickness, as far the deriving these coefficients are 
concerned, it is simply a function involved in a co-ordinate transformation. These expressions are therefore valid for 
any barotropic equilibria. 

NUMERICAL ROUTE TO A MATRIX EQUATION FOR W 

In 5j3]we arrived at the differential equation UW — by first deriving an equation for W then changed the dependent 
variable to W . Instead, we can first make the substitution W — pW and Q = pQ in Eq. [TBI — fl7l to obtain the governing 
equations for (W, Q): 



d 2 W 



d 2 W 
dZdR 



d 2 W 



dW 



Ci—- T + D! — 



dZ 2 



dR 



dW 
E 1 —+F 1 W 
oZ 



dW dW 
° 2 ~dR +E2 ~dZ + F2W + F2Q = °" 



(Bl) 
(B2) 



with 



Ar 

Fi - 

D 2 

F 2 ■■ 



ai, 



B 1 =b 1 , Ci = ci, Di=2ai 



Po 



di, Ex = bi 



PO 



2a- 



ei, 



ai- 



P'o 9 



Pa Po 9 

-- di, Ex = ex, 

-- c? 2 , E 2 = e 2 , 
h 



Cx- 



Fx=d 



P'o 
. — 
Po 



Po 
Po 

,£& 
' Po 



3 , f 
ei \-Jx, 

9 

ex h fx , 

9 



F 2 =da^ + e 2 9 - + h, 



(B3) 



We recall the unperturbed density is p = po{R)g(Z) and primes denote differentiation with respect to the argument. 
These transformation formulae make no reference to a polytropic background, so they are valid for any equilibrium 
density field separable in the above form, such as an exponential atmosphere. 
When discretized, these equations have the matrix representation 



Uxw + Uxq = 0, 
U 2 w + U 2 q = 0, 



(B4) 
(B5) 
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where q is the vector of pseudo-spectral coefficients for Q, i.e. Qi(Z) = Q(Ri, Z) = Y^,k=i Iki^kiZ / Z s ) . The matrix 
representation of UW = is then 

[U 1 -U 1 {U^ l U 2 )]w = Uw = Q. (B6) 

Note that we can divide Eq. IB2l bv F2 before the converting the operators to matrices. Then U<x is a block diagonal 
matrix of consisting only of the Chebyshev polynomials evaluated at vertical grid points. Its inverse can be pre- 
computed and stored. 

In this approach, the user only needs to specify the PDE coefficie nts defined in Appendix [5] The transformed 
coefficients A\ — F\ are used to construct the matrix U\ as described in lLinl (|2012d ). and similarly for U\ and E/jj. The 
final operator, U, results from matrix multiplication and addition, for which standard software can perform. 

ESTIMATING INSTANTANEOUS MODE GROWTH RATES 

When dealing with hydrodynamic simulations it may be impractical to frequently output data for explicit compu- 
tation of time derivatives. This is particular the case if high spatial resolution simulations are performed. However, 
we can take advantage of this and exchange time derivatives for spatial derivatives using the fluid equations. 

As usual, denote the Fourier transform with subscript m, so that 

p m (r,9,t) = / p(r,9,(j),t)exp(—im(j))d(j), (CI) 
Jo 

where we have adopted spherical co-ordinates, so here r is the spherical radius. Taking a time derivative and using 
the continuity equation gives 

^ = - / V ■ (pv) exp (-im<£)#. (C2) 







Writing this out in full, applying the usual rule for Fourier transforms to the azimuthal contribution to the divergence, 
we obtain 

dp m l<9 r2 . .-, 1 d r . . . . im 

ot r A or L J r sm 9 o9 r sin 

We can therefore just use the Fourier transform of momentum densities to calculate time derivatives of a Fourier mode. 
The complex frequency tr is defined through dtp m — iap m , from which we extract the mode frequency 10 and growth 
rate v. These are spatially-dependent when obtained from simulation data using the above procedure. So we average 
uj and v over the 9 domain and around co-rotation r S [0.8, 1.2]ro. This gives an estimate of the instantaneous growth 
rate and pattern speed of a mode with azimuthal wavenumber m at time t. 
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